!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! This module provides a unified interface to the specific test case
!! modules for the compressible Navier-Stokes code. See also \c
!! mod_testcases for implementation details.
!!
!! \n
!!
!! To add a new test case:
!! <ol>
!!
!!  <li> write a new module \c mod_testname whose public interface is
!!   given by:
!!   <ul>
!!    <li>\c mod_testname_constructor:  constructor
!!    <li>\c mod_testname_destructor:   destructor
!!    <li>\c test_name:     name of the test case (\c character)
!!    <li>\c test_description: short description (\c character)
!!    <li>\c phc: physical parameters (of type \c t_phc, refer to \c
!!    mod_physical_constants)
!!    <li>\c ntrcs: number of tracers (\c integer)
!!   </ul>
!!   Many of the following fields are optional: leave the
!!   corresponding pointers to <tt>null()</tt> if they are not
!!   implemented.
!!   <ul>
!!    <li>\c bar_mass_flow: reference mass flow
!!    \f$\underline{\mathcal{M}}_{\mathcal{I}}\f$ as discussed in the
!!    documentation of \c t_dgcomp_stv; this variable is optional,
!!    with default value zero
!!    <li>\c ref_prof:      reference profile
!!    <li>\c t_ref:         reference profile temperature (optional)
!!    <li>\c coeff_dir:     Dirichlet datum (optional, necessary for
!!    Dirichlet or open boundary conditions)
!!    <li>\c coeff_neu:     Neumann datum (optional, necessary for
!!    Neumann boundary conditions)
!!    <li>\c coeff_norm:    normal to the boundary (optional, see \c
!!    mod_dgcomp_rhs for when and how this is used)
!!    <li>\c coeff_visc:    viscosity (optional, necessary only for
!!    viscous flows)
!!    <li>\c coeff_init:    initial state
!!    <li>\c coeff_outref:  output reference
!!    <li>\c coeff_sponge:  sponge layer coefficients (optional,
!!    necessary for open boundary conditions). The output argument \c
!!    tmask can be set to <tt>.false.</tt> to indicate that none of
!!    the points in \c x requires any sponge; typically this
!!    subroutine is called once for each element, so that one can mask
!!    out some elements. There is now way to mask a single point,
!!    since this would not be very practical anayway.
!!    <li>\c coeff_f:       forcing term per unit mass (i.e.
!!    dimensionally homogeneous to an acceleration) in the momentum
!!    equation (optional)
!!   </ul>
!! 
!!  <li> use the new module in the present one;
!!
!!  <li> include the proper reference in the \c casename case blocks
!!
!!  <li> add the new module into the Makefile:
!!   <ul>
!!    <li> add the new module in the variable \c OBJ_DGTEST
!!    <li> add a target for the new module (if necessary).
!!   </ul>
!!
!! </ol>
!!
!! \todo It could be useful to introduce an additional module
!! collecting all the variables used by the test case functions, such
!! as c mass_flow in c coeff_f.
!<----------------------------------------------------------------------
module mod_dgcomp_testcases

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

 use mod_warmbubble1_test, only: &
   mod_warmbubble1_test_constructor,       &
   mod_warmbubble1_test_destructor,        &
   warmbubble1_test_name    => test_name,  &
   warmbubble1_test_description => test_description, &
   warmbubble1_phc          => phc,        &
   warmbubble1_ntrcs        => ntrcs,      &
   warmbubble1_ref_prof     => ref_prof,   &
   warmbubble1_t_ref        => t_ref,      &
   warmbubble1_coeff_init   => coeff_init, &
   warmbubble1_coeff_outref => coeff_outref

 use mod_breaking_wave_test, only: &
   mod_breaking_wave_test_constructor,     &
   mod_breaking_wave_test_destructor,      &
   breaking_wave_test_name  => test_name,  &
   breaking_wave_test_description => test_description, &
   breaking_wave_phc        => phc,        &
   breaking_wave_ntrcs      => ntrcs,      &
   breaking_wave_ref_prof   => ref_prof,   &
   breaking_wave_t_ref      => t_ref,      &
   breaking_wave_coeff_dir  => coeff_dir,  &
   breaking_wave_coeff_neu  => coeff_neu,  &
   breaking_wave_coeff_norm => coeff_norm, &
   breaking_wave_coeff_sponge => coeff_sponge, &
   breaking_wave_coeff_visc => coeff_visc, &
   breaking_wave_coeff_init => coeff_init, &
   breaking_wave_coeff_outref => coeff_outref

 use mod_boulder_windstorm_test, only: &
   mod_boulder_windstorm_test_constructor,     &
   mod_boulder_windstorm_test_destructor,      &
   boulder_windstorm_test_name  => test_name,  &
   boulder_windstorm_test_description => test_description, &
   boulder_windstorm_phc        => phc,        &
   boulder_windstorm_ntrcs      => ntrcs,      &
   boulder_windstorm_ref_prof   => ref_prof,   &
   boulder_windstorm_t_ref      => t_ref,      &
   boulder_windstorm_coeff_dir  => coeff_dir,  &
   boulder_windstorm_coeff_neu  => coeff_neu,  &
   boulder_windstorm_coeff_norm => coeff_norm, &
   boulder_windstorm_coeff_sponge => coeff_sponge, &
   boulder_windstorm_coeff_visc => coeff_visc, &
   boulder_windstorm_coeff_init => coeff_init, &
   boulder_windstorm_coeff_outref => coeff_outref

 use mod_channel_test, only: &
   mod_channel_test_constructor,     &
   mod_channel_test_destructor,      &
   channel_test_name  => test_name,  &
   channel_test_description => test_description, &
   channel_phc        => phc,        &
   channel_ntrcs      => ntrcs,      &
   channel_ref_prof   => ref_prof,   &
   channel_coeff_dir  => coeff_dir,  &
   channel_coeff_norm => coeff_norm, &
   channel_coeff_sponge => coeff_sponge, &
   channel_coeff_init => coeff_init, &
   channel_coeff_outref => coeff_outref

 use mod_turb_channel_test, only: &
   mod_turb_channel_test_constructor,     &
   mod_turb_channel_test_destructor,      &
   turb_channel_test_name  => test_name,  &
   turb_channel_test_description => test_description, &
   turb_channel_phc        => phc,        &
   turb_channel_ntrcs      => ntrcs,      &
   turb_channel_bar_mass_flow => bar_mass_flow, &
   turb_channel_ref_prof   => ref_prof,   &
   turb_channel_coeff_dir  => coeff_dir,  &
   turb_channel_coeff_visc => coeff_visc, &
   turb_channel_coeff_init => coeff_init, &
   turb_channel_coeff_outref => coeff_outref, &
   turb_channel_coeff_f => coeff_f

 use mod_straka1993_test, only: &
   mod_straka1993_test_constructor,     &
   mod_straka1993_test_destructor,      &
   straka1993_test_name  => test_name,  &
   straka1993_test_description => test_description, &
   straka1993_phc        => phc,        &
   straka1993_ntrcs      => ntrcs,      &
   straka1993_ref_prof   => ref_prof,   &
   straka1993_t_ref      => t_ref,      &
   straka1993_coeff_neu  => coeff_neu,  &
   straka1993_coeff_visc => coeff_visc, &
   straka1993_coeff_init => coeff_init, &
   straka1993_coeff_outref => coeff_outref

 use mod_gravity_wave_test, only: &
   mod_gravity_wave_test_constructor,     &
   mod_gravity_wave_test_destructor,      &
   gravity_wave_test_name  => test_name,  &
   gravity_wave_test_description => test_description, &
   gravity_wave_phc        => phc,        &
   gravity_wave_ntrcs      => ntrcs,      &
   gravity_wave_ref_prof   => ref_prof,   &
   gravity_wave_t_ref      => t_ref,      &
   gravity_wave_coeff_neu  => coeff_neu,  &
   gravity_wave_coeff_visc => coeff_visc, &
   gravity_wave_coeff_init => coeff_init, &
   gravity_wave_coeff_outref => coeff_outref

 use mod_volcano_test, only: &
   mod_volcano_test_constructor,     &
   mod_volcano_test_destructor,      &
   volcano_test_name  => test_name,  &
   volcano_test_description => test_description, &
   volcano_phc        => phc,        &
   volcano_ntrcs      => ntrcs,      &
   volcano_ref_prof   => ref_prof,   &
   volcano_t_ref      => t_ref,      &
   volcano_coeff_dir  => coeff_dir,  &
   volcano_coeff_neu  => coeff_neu,  &
   !volcano_coeff_sponge => coeff_sponge, &
   volcano_coeff_visc => coeff_visc, &
   volcano_coeff_init => coeff_init, &
   volcano_coeff_outref => coeff_outref

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dgcomp_testcases_constructor, &
   mod_dgcomp_testcases_destructor,  &
   mod_dgcomp_testcases_initialized, &
   test_name,   & ! name of the test case
   test_description,& ! short description of the test case
   t_phc, phc,      & ! physical constants
   ntrcs,           & ! number of tracers
   bar_mass_flow,   & ! reference mass flow
   ref_prof, t_ref, & ! reference profile
   coeff_dir,   & ! Dirichlet datum
   coeff_neu,   & ! Neumann datum
   coeff_norm,  & ! boundary normal
   coeff_visc,  & ! viscosity
   coeff_init,  & ! initial state for time dependent problems
   coeff_outref,& ! output reference
   coeff_sponge,& ! sponge coefficients
   coeff_f        ! forcing term

 ! Compiler bug: this is DPD200244054
 !  http://software.intel.com/en-us/forums/topic/392251
 ! The private blanket attribute should be restored as soon as
 ! possible.
 !private

!-----------------------------------------------------------------------

 ! Module types and parameters
 integer, protected :: ntrcs ! used in the interfaces
 abstract interface
  pure function i_coeff_dir(x,u_r,u_i,breg) result(ub)
   import :: wp, ntrcs
   implicit none
   real(wp), intent(in) :: x(:,:)
   real(wp), intent(in), optional :: u_r(:,:), u_i(:,:)
   integer, intent(in), optional :: breg
   real(wp) :: ub(2+size(x,1)+ntrcs,size(x,2))
  end function i_coeff_dir
 end interface
 abstract interface
  pure function i_coeff_neu(x,u_r,breg) result(h)
   import :: wp, ntrcs
   implicit none
   real(wp), intent(in) :: x(:,:)
   real(wp), intent(in), optional :: u_r(:,:)
   integer, intent(in), optional :: breg
   real(wp) :: h(size(x,1),1+size(x,1)+ntrcs,size(x,2))
  end function i_coeff_neu
 end interface
 abstract interface
  pure function i_coeff_norm(x,breg) result(nb)
   import :: wp
   implicit none
   real(wp), intent(in) :: x(:,:)
   integer, intent(in) :: breg
   real(wp) :: nb(size(x,1),size(x,2))
  end function i_coeff_norm
 end interface
 abstract interface
  pure function i_coeff_visc(x,rho,p,t) result(nu)
   import :: wp
   implicit none
   real(wp), intent(in) :: x(:,:), rho(:), p(:), t(:)
   real(wp) :: nu(2,size(x,2))
  end function i_coeff_visc
 end interface
 abstract interface
  pure function i_coeff_init(x,u_r) result(u0)
   import :: wp, ntrcs
   implicit none
   real(wp), intent(in) :: x(:,:), u_r(:,:)
   real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))
  end function i_coeff_init
 end interface
 abstract interface
  pure subroutine i_coeff_sponge(tmask,tau,taup,x)
   import :: wp
   implicit none
   real(wp), intent(in) :: x(:,:)
   logical, intent(out) :: tmask
   real(wp), intent(out) :: tau(:), taup(:)
  end subroutine i_coeff_sponge
 end interface
 abstract interface
  pure function i_coeff_f(x,mass_flow,imass_flow) result(f)
   import :: wp
   implicit none
   real(wp), intent(in)  :: x(:,:)
   real(wp), intent(in) :: mass_flow(:), imass_flow(:)
   real(wp) :: f(size(x,1),size(x,2))
  end function i_coeff_f
 end interface

 ! Module variables
 ! public members
 character(len=10000), protected              :: test_name
 character(len=10000), allocatable, protected :: test_description(:)
 type(t_phc), save, protected                 :: phc
 real(wp), allocatable, protected             :: bar_mass_flow(:)
 character(len=10000), protected              :: ref_prof
 real(wp), pointer, protected                 :: t_ref => null()
 procedure(i_coeff_dir   ), pointer :: coeff_dir    => null()
 procedure(i_coeff_neu   ), pointer :: coeff_neu    => null()
 procedure(i_coeff_norm  ), pointer :: coeff_norm   => null()
 procedure(i_coeff_visc  ), pointer :: coeff_visc   => null()
 procedure(i_coeff_init  ), pointer :: coeff_init   => null()
 procedure(i_coeff_init  ), pointer :: coeff_outref => null()
 procedure(i_coeff_sponge), pointer :: coeff_sponge => null()
 procedure(i_coeff_f     ), pointer :: coeff_f      => null()
 protected :: coeff_dir, coeff_neu, coeff_norm, coeff_visc, &
   coeff_init, coeff_outref, coeff_sponge, coeff_f     

 logical, protected ::               &
   mod_dgcomp_testcases_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_dgcomp_testcases'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_dgcomp_testcases_constructor(testname,d)
  character(len=*), intent(in) :: testname
  integer, intent(in) :: d

  integer :: i
  character(len=10000) :: message
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized          .eqv..false.) .or. &
       (mod_kinds_initialized             .eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dgcomp_testcases_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! set the defaults
   allocate(bar_mass_flow(d)); bar_mass_flow = 0.0_wp

   casename: select case(testname)

    case(warmbubble1_test_name)
     call mod_warmbubble1_test_constructor(d)
     test_name = warmbubble1_test_name
     allocate(test_description(size(warmbubble1_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(warmbubble1_test_description(i))
     enddo
     phc   = warmbubble1_phc
     ntrcs = warmbubble1_ntrcs
     ref_prof = trim(warmbubble1_ref_prof)
     t_ref        => warmbubble1_t_ref
     coeff_init   => warmbubble1_coeff_init
     coeff_outref => warmbubble1_coeff_outref

    case(breaking_wave_test_name)
     call mod_breaking_wave_test_constructor(d)
     test_name = breaking_wave_test_name
     allocate(test_description(size(breaking_wave_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(breaking_wave_test_description(i))
     enddo
     phc   = breaking_wave_phc
     ntrcs = breaking_wave_ntrcs
     ref_prof = trim(breaking_wave_ref_prof)
     t_ref        => breaking_wave_t_ref
     coeff_visc   => breaking_wave_coeff_visc
     coeff_init   => breaking_wave_coeff_init
     coeff_outref => breaking_wave_coeff_outref
     coeff_dir    => breaking_wave_coeff_dir
     coeff_neu    => breaking_wave_coeff_neu
     coeff_norm   => breaking_wave_coeff_norm
     coeff_sponge => breaking_wave_coeff_sponge

    case(boulder_windstorm_test_name)
     call mod_boulder_windstorm_test_constructor(d)
     test_name = boulder_windstorm_test_name
    allocate(test_description(size(boulder_windstorm_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(boulder_windstorm_test_description(i))
     enddo
     phc   = boulder_windstorm_phc
     ntrcs = boulder_windstorm_ntrcs
     ref_prof = trim(boulder_windstorm_ref_prof)
     t_ref        => boulder_windstorm_t_ref
     coeff_visc   => boulder_windstorm_coeff_visc
     coeff_init   => boulder_windstorm_coeff_init
     coeff_outref => boulder_windstorm_coeff_outref
     coeff_dir    => boulder_windstorm_coeff_dir
     coeff_neu    => boulder_windstorm_coeff_neu
     coeff_norm   => boulder_windstorm_coeff_norm
     coeff_sponge => boulder_windstorm_coeff_sponge

    case(channel_test_name)
     call mod_channel_test_constructor(d)
     test_name = channel_test_name
     allocate(test_description(size(channel_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(channel_test_description(i))
     enddo
     phc   = channel_phc
     ntrcs = channel_ntrcs
     ref_prof = trim(channel_ref_prof)
     coeff_init   => channel_coeff_init
     coeff_outref => channel_coeff_outref
     coeff_dir    => channel_coeff_dir
     coeff_norm   => channel_coeff_norm
     coeff_sponge => channel_coeff_sponge

    case(turb_channel_test_name)
     call mod_turb_channel_test_constructor(d)
     test_name = turb_channel_test_name
     allocate(test_description(size(turb_channel_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(turb_channel_test_description(i))
     enddo
     phc   = turb_channel_phc
     ntrcs = turb_channel_ntrcs
     bar_mass_flow = turb_channel_bar_mass_flow
     ref_prof = trim(turb_channel_ref_prof)
     coeff_visc   => turb_channel_coeff_visc
     coeff_init   => turb_channel_coeff_init
     coeff_outref => turb_channel_coeff_outref
     coeff_dir    => turb_channel_coeff_dir
     coeff_f      => turb_channel_coeff_f

    case(straka1993_test_name)
     call mod_straka1993_test_constructor(d)
     test_name = straka1993_test_name
     allocate(test_description(size(straka1993_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(straka1993_test_description(i))
     enddo
     phc   = straka1993_phc
     ntrcs = straka1993_ntrcs
     ref_prof = trim(straka1993_ref_prof)
     t_ref        => straka1993_t_ref
     coeff_visc   => straka1993_coeff_visc
     coeff_init   => straka1993_coeff_init
     coeff_outref => straka1993_coeff_outref
     coeff_neu    => straka1993_coeff_neu

    case(gravity_wave_test_name)
     call mod_gravity_wave_test_constructor(d)
     test_name = gravity_wave_test_name
     allocate(test_description(size(gravity_wave_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(gravity_wave_test_description(i))
     enddo
     phc   = gravity_wave_phc
     ntrcs = gravity_wave_ntrcs
     ref_prof = trim(gravity_wave_ref_prof)
     t_ref        => gravity_wave_t_ref
     coeff_visc   => gravity_wave_coeff_visc
     coeff_init   => gravity_wave_coeff_init
     coeff_outref => gravity_wave_coeff_outref
     coeff_neu    => gravity_wave_coeff_neu

    case(volcano_test_name)
     call mod_volcano_test_constructor(d)
     test_name = volcano_test_name
    allocate(test_description(size(volcano_test_description)))
     do i=1,size(test_description)
       test_description(i) = trim(volcano_test_description(i))
     enddo
     phc   = volcano_phc
     ntrcs = volcano_ntrcs
     ref_prof = trim(volcano_ref_prof)
     t_ref        => volcano_t_ref
     coeff_visc   => volcano_coeff_visc
     coeff_init   => volcano_coeff_init
     coeff_outref => volcano_coeff_outref
     coeff_dir    => volcano_coeff_dir
     coeff_neu    => volcano_coeff_neu
     !coeff_sponge => volcano_coeff_sponge

    case default
     write(message,'(a,a,a)') 'Unknown test case "',trim(testname),'"'
     call error(this_sub_name,this_mod_name,message)
   end select casename

   write(message,'(a,a,a)') 'Test case: "',trim(test_name),'"'
   call info(this_sub_name,this_mod_name,message)

   mod_dgcomp_testcases_initialized = .true.
 end subroutine mod_dgcomp_testcases_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dgcomp_testcases_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dgcomp_testcases_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   casename: select case(trim(test_name))
    case(warmbubble1_test_name)
     call mod_warmbubble1_test_destructor()
    case(breaking_wave_test_name)
     call mod_breaking_wave_test_destructor()
    case(boulder_windstorm_test_name)
     call mod_boulder_windstorm_test_destructor()
    case(channel_test_name)
     call mod_channel_test_destructor()
    case(turb_channel_test_name)
     call mod_turb_channel_test_destructor()
    case(straka1993_test_name)
     call mod_straka1993_test_destructor()
    case(gravity_wave_test_name)
     call mod_gravity_wave_test_destructor()
    case(volcano_test_name)
     call mod_volcano_test_destructor()
   end select casename

   ! nullify here the procedure pointers
   nullify( coeff_dir,coeff_neu,coeff_norm,coeff_visc,coeff_init, &
     coeff_outref,coeff_sponge,coeff_f )

   deallocate(test_description,bar_mass_flow)

   mod_dgcomp_testcases_initialized = .false.
 end subroutine mod_dgcomp_testcases_destructor

!-----------------------------------------------------------------------

end module mod_dgcomp_testcases

